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Circadian (~24h) timekeeping is essential for the lives of many organisms. To understand the 
biochemical mechanisms of this timekeeping, we have developed a detailed mathematical model of 
the mammalian circadian clock. Our model can accurately predict diverse experimental data 
including the phenotypes of mutations or knockdown of clock genes as well as the time courses and 
relative expression of clock transcripts and proteins. Using this model, we show how a universal 
motif of circadian timekeeping, where repressors tightly bind activators rather than directly binding 
to DNA, can generate oscillations when activators and repressors are in stoichiometric balance. 
Furthermore, we find that an additional slow negative feedback loop preserves this stoichiometric 
balance and maintains timekeeping with a fixed period. The role of this mechanism in generating 
robust rhythms is validated by analysis of a simple and general model and a previous model of the 
Drosophila circadian clock. We propose a double-negative feedback loop design for biological clocks 
whose period needs to be tightly regulated even with large changes in gene dosage. 
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Introduction 

Orcadian (~24h) clocks time many physiological and meta- 
bolic processes. When these clocks were first discovered, 
three basic properties were identified (Dunlap et al 2004). 
(1) Rhythms need to be autonomous. (2) Rhythms need to 
be capable of adjusting in response to external signals. 
(3) Rhythms need to persist over a wide range of temperatures. 
More recently, the biochemical mechanisms of circadian 
timekeeping have been identified (Ko and Takahashi, 2006). 
In particular, interlocked transcription-translation feedback 
loops (TTFLs) have been discovered as the basic mechanism of 
rhythm generation in many organisms (Novak and Tyson, 
2008) . With this discovery, recent experimentation has identi- 
fied another property of circadian rhythms in higher organisms. 
Circadian rhythms persist with a 24-h period even in the 
presence of large changes in the expression of the components 
of these TTFLs (Ko and Takahashi, 2006; Dibner et al 2009). 
While mechanisms for rhythm generation with a flexible period 
have been identified (Strieker et al 2008; Tsai et al, 2008; Tigges 
et al 2009), mechanisms for this robustness of period to gene 
dosage remain unexplained, even by mathematical models 
(Dibner et al 2009). 

Two interlocked negative feedback loops have been identi- 
fied in the TTFL networks generating circadian rhythms in 
higher organisms (Figure 1) (Blau and Young, 1999; Glossop 
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etal 1999; Benito etal 2007; Liu et al 2008). A 'core' negative 
feedback loop consists of repressors (PERIOD and TIMELESS 
in Drosophila or PERIOD1-3 and CRYPTOCHROME1-2 in 
mammals), which inactivate activators (CYCLE and CLOCK in 
Drosophila and BMAL1-2 and CLOCK in mammals) of their 
own transcription. An additional negative feedback loop 
controls the expression of the activators, which inactivate 
their own transcription through Vrille {Drosophila) or the Rev- 
erbs genes (Mammals) (Blau and Young, 1999; Preitner et al 
2002). While other feedback loops have also been identified, 
these two negative feedback loops seem to predominate (Blau 
and Young, 1999; Glossop et al 1999; Benito et al 2007; Liu 
et al 2008; Bugge et al 2012; Cho et al 2012). 

Near 24-h oscillations persist even when the components of 
the TTFLs of the circadian clock are over or under expressed. 
Heterozygous mutations of clock genes never abolish rhyth- 
micity, and their period phenotypes are either indistinguish- 
able from the wild-type (WT) phenotypes or much smaller 
than mutations that affect post-translational modifications 
(Baggs et al 2009; Etchegaray et al 2009; Lee et al 2009). 
Abolishing rhythmicity through single gene knockout is 
surprisingly difficult (Baggs et al 2009; Ko et al 2010). 
Moreover, the mammalian circadian clock is also resistant to 
global changes in transcription rates (Dibner et al 2009). 
These results all suggest that gene dosage may not be 
important for circadian timekeeping in higher organisms. 
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Gene dosage, however, is not completely unimportant 
for timekeeping. Knockdown of clock genes causes 
increased expression in similar components through paralog 
compensation, which may help restore gene dosage and 
indicates that gene dosage needs to be tightly regulated (Baggs 
et al, 2009). Population rhythmicity in mouse embryonic 
fibroblasts shows much lower amplitude than in liver, which 
might be due to the fact that the ratio of repressors to activators 
is significantly lower in fibroblasts than that found in liver (Lee 
etal, 2001, 2011). A 1-1 stoichiometric binding occurs between 
the activators and repressors driving rhythms in Drosophila 
(Menet et al, 2010), although not in Neurospora (He et al, 2005; 
Huang et al, 2007). 

Here, we propose a mechanistic explanation for the 
robustness to gene dosage in the circadian clock of higher 
organisms through mathematical modeling. We develop the 
most detailed mathematical model of the mammalian circa- 
dian clock available, which should be useful in many future 
studies. Our model reproduces a surprising amount of 
experimental data on the mammalian circadian clock includ- 
ing the time courses and relative concentrations of key 
transcripts and proteins, the effects of mutations of key clock 
genes, and the effects of changes in gene dosage. With this 
model, we show that proper stoichiometric balance between 
activators (BMAL-CLOCK/NPAS2) and repressors (PERI -2/ 
CRY1-2) is key to sustained oscillations. Furthermore, we find 
that an additional slow negative feedback loop, in which 
activators indirectly inactivate themselves, improves the 
regulation of the stoichiometric balance and sustains oscilla- 
tions with a nearly constant period over a large change in gene 
expression level. Tight binding between activators and 
repressors is also predicted to be crucial for rhythm genera- 
tion. These mechanisms are also validated by mathematical 
analysis of a simplified mathematical model of the mammalian 
circadian clock, and simulations of a previously published 
Drosophila model. We here propose a novel design for 
biological oscillators where maintaining period is crucial: a 
core negative feedback loop with repression by protein 
sequestration, with an additional negative feedback loop, 
which controls a relatively stable activator. 

Results 

Mathematical modeling of the mammalian 
circadian clock 

We develop a new mathematical model of the intracellular 
mammalian circadian clock. This model contains key genes, 
mRNAs and proteins (PERI, PER2, CRY1, CRY2, BMAL1/2, 
NPAS2, CLOCK, CKIe/5, GSK3(3, Rev-erba/(3) that have been 
found to be central to mammalian circadian timekeeping 
(Figure 1A). While greatly expanded, the model is largely 
based on our previous model, which has made surprising 
predictions about mammalian timekeeping that have been 
subsequently verified experimentally (Forger and Peskin, 
2003; Gallego et al, 2006; Ko et al, 2010; Yamada and Forger, 
2010) . Modifications and extensions of the model are described 
in the Materials and methods, Supplementary information and 
Supplementary Tables 1 and 2. The parameters of the model 
are estimated using experimental data and a simulated 
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annealing method (a global stochastic parameter searcher) 
(Gonzalez et al, 2007) (see Materials and methods, 
Supplementary information and Supplementary Table 3 for 
details). In particular, we incorporated experimentally deter- 
mined rate constants (Supplementary Table 3) (Kwon et al, 
2006; Siepka et al, 2007; Chen et al, 2009; Suter et al, 2011), fit 
the time courses of both mRNA and proteins (Figure 2A and B) 
(Lee et al, 2001; Reppert and Weaver, 2001; Ueda et al, 2005) 
and fit the relative abundance of proteins (Figure 2C) (Lee 
et al, 2001). 

Our model accurately predicts the phenotype of known 
mutations of genes in the central circadian clock (suprachias- 
matic nuclei, SCN) (Yoo et al, 2005; Baggs et al, 2009; Ko et al, 
2010), which other models do not predict (Table I) (Forger and 
Peskin, 2003; Leloup and Goldbeter, 2003; Mirsky et al, 2009; 
Relogio et al, 2011). Interestingly, our model shows opposite 
phenotypes for Cryl ~ f ~ and Cry 2 ~ 1 ~ matching experimental 
data (Liu et al, 2007) . There are two differences between CRY1 
and CRY2 in our model. First, Cryl transcription is delayed 
through repression by Rev-erba and Rev-erb(3 (Preitner et al, 
2002; Liu etal, 2008; Ukai-Tadenuma etal, 2011). Additionally, 
Cryl mRNA is more stable than Cry2 mRNA and CRY1 protein 
is more stable than CRY2 protein (Busino et al, 2007; Siepka 
et al, 2007; Chen et al, 2009). Since a longer half-life causes 
rhythms to be delayed, and delayed rhythms cause a longer 
period (Forger, 2011; Ukai-Tadenuma et al, 2011), removing 
CRY1 shortens the period and removing CRY2 lengthens the 
period. The opposite phenotypes of Clock~ f ~ (null mutation) 
and Clock A19/+ (dominant-negative mutation) are also cor- 
rectly simulated in the model for the first time (Vitaterna et al, 
1994; Herzog etal, 1998; Debruyne etal, 2006). Moreover, our 
model also predicts the mutant phenotypes of isolated SCN 
neurons, which are different from the SCN slices (Liu et al, 
2007). We note that SCN slices have significantly higher gene 
expression of perl and per2 through CREB/CRE pathway than 
isolated SCN neurons (Yamaguchi et al, 2003). Interestingly, 
when we reduced perl and per2 expression about 60% in our 
model, our model was able to accurately reproduce the 
phenotypes of isolated SCN neurons (Table II) . 

We also conducted a sensitivity analysis to look at what 
parameters determine the period of our model. Four of the top 
five high parameters, in our sensitivity analysis, were also in 
the top five found in a previous sensitivity analysis with the 
original Forger and Peskin model and which was used to 
conclude that PER2 plays a dominant role in period 
determination (Wilkins et al, 2007) (see Supplementary 
Figure 1). 



Proper stoichiometric balance between activators 
and repressors is crucial to sustained rhythms 

Since our mathematical model can accurately predict the 
phenotype of known mutations of the mammalian circadian 
clock, we next looked for a mechanism that could explain why 
some phenotypes were rhythmic, while others were not. We 
found that stoichiometry plays a key role in determining which 
mutations showed rhythmic phenotypes. Here, we define 
stoichiometry as the average ratio between the concentration 
of repressors (all forms of PER and CRY in the nucleus) to that 
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Figure 1 Schematic of the detailed mammalian circadian clock model. (A) Only some of the relevant species are shown. Circles refer to transcripts and squares are 
proteins, possibly in complex. Small circles refer to phosphorylation states that are color coded by the kinases that perform the phosphorylation. See section 'Description 
of the detailed model' in Supplementary information for details. (B) The detailed model consists of a core negative feedback loop and an additional negative feedback 
loop (the NNF structure). The repressors (PER1-2 and CRY1-2) inactivate the activators (BMALs and CLOCK/NPAS2) of their own transcription expression through the 
core negative feedback loop. The activators inactivate their own transcription expression by inducing the Rev-erbs through the secondary negative feedback loop. 



of activators (all forms of BMAL-CLOCK/NPAS2 in the 
nucleus) over a period. Moreover, we specifically refer to 
repressors and activators of E/E'-boxes when discussing 
stoichiometry. We found that mutations that caused the 
stoichiometry to be too high or too low, yielded arrhythmic 
phenotypes (Figure 3 A) . So long as the mutations allowed the 
stoichiometry to be around a 1-1 ratio, relatively high 
amplitude oscillations were seen. Thus, we predict that 
stoichiometry provides a unifying principle to determine the 
rhythmicity of mutations of the mammalian circadian 
clock. To further test this principle, we constitutively expressed 
either the Per2 gene (the dominant repressor gene) or the 
Bmal and Clock genes (the dominant activator genes) 
at different levels. Interestingly, within a range centered near 



a 1-1 stoichiometry, the model shows sustained oscillations 
with high amplitude (Figure 3B) . However, if the stoichiometry 
was too high or too low, rhythms are dampened or completely 
absent (Figure 3B). This matches a recent experimental study 
showing that the amplitude and sustainability of population 
rhythms increase when the level of PER-CRY is increased 
closer to that of BMAL1 -CLOCK in mouse fibroblasts (Lee 
et al, 2011). 

We defined the stoichiometry as the average ratio between 
the total concentrations of repressors to that of activators 
over a period. However, recent work has shown that CRY1 
has stronger repressor activity than CRY2. The underlying 
biochemical mechanisms for this result have not been fully 
identified (Khan et al, 2012). If the difference is due to a 
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Figure 2 Validation of the detailed model. (A) Predicted mRNAs time courses in SCN (Ueda et al, 2005). Time courses were normalized so that the peak value is 1 , 
matching experimental data. (B) Predicted protein time courses in SCN (Reppert and Weaver, 2001). As had been done previously, we normalize the protein time 
courses so that the maximum is 1 and the minimum is 0. (C) Model comparison of the relative abundance of proteins in liver and fibroblast (Lee et al, 2001 , 2009, 201 1 ). 
All of the values were normalized so that the maximum abundance of the CRY1 protein is 1. For the CKIs/5, CKIs maximal expression is -22.5% of the maximum 
abundance of CRY1 in the liver (Lee et al, 2001) and CKIS is two times more abundant than CKIs in the fibroblast (Lee etal, 2009). From this, we assumed that total 
CKIe/5 would be -67.5% of the maximum value of CRY1 in mice liver and fibroblast. 



different post-translational mechanism (e.g., binding between 
PER and CRY, which could affect the repressor concentration in 
the nucleus), the current definition of stoichiometry can be 
kept. Otherwise, a more sophisticated definition of stoichio- 
metry may be needed (e.g., one that gives more weight to 
concentration of CRY1 than that of CRY2). 



How stoichiometry generates rhythms 

To test the role of stoichiometry in sustaining oscillations, we 
developed a simple model by modifying the well-studied 
Goodwin model (Goodwin, 1965) to include an activator (A), 
which becomes inactive when bound by a repressor (P) 
(Figure 3C). Transcription is proportional to the fraction of 
free activator that is not bound by the repressor, f(P, A, K d ) 
(Buchler and Cross, 2009), matching experimental data from 
the mammalian circadian clock (Supplementary Figure 2) 
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(Froy et al, 2002). mRNA (M) is translated to a repressor 
protein (P c ) . The protein enters the nucleus (P) and binds and 
inhibits the activator (A). This generates a single-negative 
feedback loop (SNF) since the activator is constitutively 
expressed. The model is similar to a previously published 
mathematical model (Francois and Hakim, 2005); however, 
we allow for both association and dissociation of the activator 
and repressor (through a defined K d ), which turns out to be 
crucial for understanding the effects of stoichiometry. By 
nondimensionalization and setting the clearance rates of all 
species to be equal (to increase the chance of oscillations, 
see Forger, 2011), only two parameters remain: the activator 
concentration (A) and the dissociation constant {K d ) (see 
Supplementary information) . 

When we changed the activator concentration, which 
changed the stoichiometry (average ratio between the level 
of repressor (P) to the level of activator (A)), sustained 
oscillations were only seen at around a 1-1 stoichiometry 
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Table I Comparison of model predictions with experimental data and previous model predictions on the phenotypes of circadian mutations 











Relogio et ol 


Mirsky 


Leloup and 


Forger and 


Gene 




i 

Animal 


A 1 

New model 




cl CLL l^vvy ) 


LrOiQueier ^zuuo j 


i esKin ^zuuo j 


f n; 7 - / - 


oiiui i 


oiiui i 


_ ^ 


Long 


AR 




WT 




Long 


Long 


i i £ 

-f- 1.0 


Long 


Long 




Long 


Perl~'~ 






WT 


AR 


AR 


Short 


Long 


Perl ldc 


WT 


Short/ AR 










Per2~ / ~ 






AR 


AR 


AR 


Short 


Short 


Per2 ldc 




Short/ AR 












Bmall~'~ 


SR** 




AR 


AR 


AR 


AR 


AR 


Bmall~ / + 


WT* 




+ 0.1 


AR 


NA 


AR 


Long 


Clock- 


WT 


Short 


-0.2 


Long 


AR 


AR 


AR 


Clock Al9/A19 


AR* 


Long 


AR 


Long 


NA 


NA 


NA 


Clock A19/ + 


Long* 




+ 1.1 


Long 


NA 


NA 


NA 


Npas2 1 


WT 


Short 


WT 


NA 


NA 


NA 


NA 


Rev-erba ~'~ 




Short 


-0.2 


AR 


Short 


NA 


WT 


CKlE tau/tau 


Short 


Short 


-3 


NA 


NA 


Short 


Short 



Here we indicate whether the phenotype predicted by our model, or seen in experimental data is WT, stochastically rhythmic (SR), arrhythmic (AR) or shows a change 
in period in hours. Experimental data can be found in Baggs et al (2009) as well as references cited therein, except those marked with * which can be found in Yoo et al 
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new model. NA represents not available. For the Leloup-Goldbeter model, first parameter set of the model is used. 



similar to our detailed model (Figure 3D). As the other 
parameter (K d ) decreased (indicating tight binding), the range 
of stoichiometry that permitted oscillation increased 
(Figure 3D). Interestingly, if the binding was too weak, the 
rhythms did not occur. The tight binding between activators 
and repressors is also found in the detailed model, and in the 
mammalian circadian clock (Lee et al, 2001; Froy et al, 2002; 
Sato et al, 2006). This indicates that the sustained rhythms 
require tight binding as well as balanced stoichiometry in the 
circadian clock. 

Many previous studies have argued that ultrasensitive 
responses (e.g., a large change in transcription rate for a small 
change in repressor or activator concentration) can cause 
oscillations in feedback loops (Kim and Ferrell, 2007; Buchler 
and Louis, 2008; Novak and Tyson, 2008; Forger, 2011). A 
previous study showed that an ultrasensitive response can be 
generated by tight binding of activators and repressors in a 
synthetic system (Buchler and Cross, 2009). Taken together, 
this provides a potential mechanism of rhythm generation. 
That is, when the total concentration of repressor is higher 
than that of activators, the repressor sequesters and buffers 
activator and inhibits transcription completely (Buchler and 
Louis, 2008). As the repressor is depleted, the excess free 
activators are no longer sequestered by repressors and are free 
to turn on the transcription. At this threshold, transcription of 
repressor shows an ultrasensitive response to the concentra- 
tion of repressor or activator. Ultrasensitive responses amplify 
rhythms and prevent rhythms from dampening (Forger, 2011). 
In both our simple and our detailed model, we found 
ultrasensitive responses around a 1-1 stoichiometry 
(Supplementary Figure 3 A) . When the stoichiometry was not 
around 1-1, an ultrasensitive response was not seen, and both 
models did not show sustained rhythms. 

Over the course of a day, as levels of repressor and activator 
change, the stoichiometry and also sensitivity change as well. 
We found that the 1-1 average stoichiometry is required to 
generate the ultrasensitive response, which causes rhythms 
through mathematical analysis, confirming our simulation 
results (Figure 3D). That is, via both local and global stability 
analysis, we derived an approximate range of the 



Table II Comparison of modified model predictions with experimental data of 
single SCN neurons on the phenotypes of circadian mutations 



Gene 


dSCN 


Model 


Cryl-i- 


AR 


AR 


Cry2-'- 


Long 


+ 2.3 


Perl~'~ 




AR 


Perl ldc 


AR 




Bmall-l- 


AR* 


AR 



Here, we indicate whether the phenotype predicted by our model, or seen in 
experimental data is arrhythmic (AR) or shows a change in period in hours. 
Experimental data can be found in Liu et al (2007), except those marked with * 
which can be found in Ko et al (2010). See Materials and methods for details. 



stoichiometries (<S>) that permit oscillations 

g<(S)< * 

(see Supplementary information). In agreement with our 
simulations shown in Figure 3D, this mathematical analysis 
also suggests that (1) oscillations are seen around a 1-1 
stoichiometry; (2) the stoichiometry needs to be >8/9 for 
sustained rhythmicity; (3) as the binding between acti- 
vators and repressors becomes tighter, the upper bound on 
stoichiometry increases; (4) if the binding is too weak 
(e.g., K d = 10 ~ 3 ), sustained oscillations do not occur. 



An additional negative feedback loop improves the 
regulation of stoichiometric balance 

If stoichiometry is key to sustained oscillation, are there 
mechanisms within circadian clocks that keep the stoichio- 
metry of components balanced? Does the additional negative 
feedback loop of the negative-negative feedback loop (NNF) 
structure, found in circadian clocks, help balance stoichiome- 
try? To test this structure, we added an additional negative 
feedback loop into our simple model (Figure 4A). Previously, 
other studies suggested that an additional positive, rather than 
negative, feedback loop could sustain intracellular clocks 
(Barkai and Leibler, 2000; Strieker et al, 2008; Tsai et al, 2008; 
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Figure 3 Proper stoichiometry between activators and repressors is the key to sustained oscillations. (A) Our detailed mathematical model accurately predicts the 
phenotype of the known mutations in circadian genes (Table I). We plot the stoichiometry predicted by our model in these mutants with the relative amplitude of Perl 
mRNA rhythms (or Per2 mRNA when considering the Perl ~ '~ ). Here, an amplitude of zero means rhythms are not sustained. These results indicate that the phenotype 
of the mutants can be predicted by their effects on stoichiometry. (B) The stoichiometry between repressors and activators is changed by constitutively expressing either 
the Per2 gene or the Bmals and Clock genes at different levels. Note that the model is rhythmic only when the stoichiometry is near 1-1 . The relative amplitude of the 
Perl mRNA is measured. (C) Schematic of a simplified model based on the Goodwin oscillator. Instead of a Hill-type equation, the sequestration of the activator (A) by 
the repressor (P) is used to describe repression of the gene. (D) Oscillations are seen around a 1-1 stoichiometry as the level of activator is changed. The range of the 
stoichiometry widens as the dissociation constant (/C d ) decreases or the binding between the activator and the repressor tightens. 



Tigges et al, 2009). We tested these structures by including an 
additional protein R (Rev-ERBs or RORs in the mammalian 
circadian clocks) that is transcribed in a similar way to P. R 
then represses (as in the Rev-erbs) or promotes (as in the Rors) 
the production of A in the negative-negative feedback loop 
(NNF) or the positive-negative feedback loop (PNF) structure, 
respectively (Figure 4A) . 

We studied how the SNF, NNF and PNF structures effectively 
maintain the stoichiometric balance when model parameters 
(e.g., transcription rate) are changed. With both simulation 
and steady-state analysis, we found that the NNF structure is 
best at keeping stoichiometry balanced while the PNF 
structure is worst at keeping stoichiometry balanced, regard- 
less which parameters are perturbed (see Supplementary 
information, Figure 4B and Supplementary Figure 4A-C). 
Moreover, our detailed model, which also follows the NNF 
structure, also carefully balanced the stoichiometry by 
controlling the expression of repressors and activators. 
Knockdown of the repressor Cryl leads to higher expression 
of the repressors, which are controlled by E-boxes, and lower 
expression of the activators, which are controlled by a ROREs 
(Figure 4C). Opposite effects are seen when the activator 
CLOCK is removed (Figure 4C). This active control of 
repressors and/or activators via the NNF structure regulates 

6 Molecular Systems Biology 2012 



the stoichiometric balance tightly (Supplementary Figure 4D) 
and matches experimental data on gene dosage (Baggs et al, 
2009). Moreover, the detailed model (with the NNF structure) 
also correctly predicts the change of clock gene expression 
after the removal of the additional negative feedback loop 
(£ev-erboc,(3 _/ ~) (Figure4D) (Liu etal, 2008; Bugge etal, 2012; 
Cho et al, 2012). In particular, knockout of the Rev-erbai,fi 
decreases PER expression, but increase CRY1 expression. For 
our nominal set of parameters, oscillations are still possible 
when this additional negative feedback is removed. However, 
for other sets of parameters, where stoichiometry is not as well 
balanced, removal of this additional negative feedback stops 
rhythmicity (see below). This could explain the phenotype of 
the Rev-erbu,$~ / ~ , which show some indications of rhythmi- 
city (Bugge et al, 2012; Cho et al, 2012). Our model predicts 
that rhythm generation remains in cell types that have a near 
balanced stoichiometry, and a lack of rhythms in cell types 
without a balanced stoichiometry. 

A slow additional negative feedback loop 
improves the robustness of rhythms 

Our central hypothesis is that, as stoichiometry is more tightly 
regulated, oscillations will occur over a wider range of 
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Figure 4 The NNF structure maintains stoichiometry in balance by active compensation of both repressors and activators. (A) A negative or positive feedback 
controlling the activator is added to the original negative feedback controlling the repressor. (B) The relative sensitivity (% change in mean level of stoichiometry per % 
change in transcription rate of repressor) in the simple models with SNF, NNF and PNF structure were measured over a range of the transcription rates of repressor (see 
Supplementary Figure 4A). Then, we calculated the average of relative sensitivity over the range of parameters. On average, the relative sensitivity of the NNF model is 
about two-fold less sensitive than that of the SNF model, but that of the PNF model is about four-fold more sensitive than that of the SNF model (see Supplementary 
information and Supplementary Figure 4A-C for details). (C) The detailed model matches data from Gene Dose Network Analysis experiments (Baggs era/, 2009). After 
the knockout of a repressor gene (here, Cry1), the activity of the repressor promoters, controlled by an E-box, increases. This increases the expression of Rev-Erbs and 
reduces the activity of the activator promoter, controlled by a RORE. An opposite phenotype is seen when an activator (here, Clock) is knocked out. The activity of E-box 
decreases. This decreases the expression of Rev-Erbs and increases the activity of RORE. This active compensation through the NNF structure allows the stoichiometry 
to be balanced after the repressors or activators knockout. (D) The detailed model matches data from Rev-erbs''' (Liu era/, 2008; Bugge era/, 2012; Cho era/, 2012). 
Rev-erba''' (50% reduction of transcription rate of the Rev-erbs due to the presence of Rev-erbp) slightly shortens the period and has little effect on the expression 
level of Per2, Cry1 and Bmall. Double knockout of the Rev-erba and Rev-erbp (100% reduction of transcription rate of the Rev-erbs) increases the expression level of 
Bmall and Cry1, but decreases that of Per2. All the values were normalized by the average of Per2 expression level in WT. 



parameters. To confirm this, we varied the transcription rate of 
the activator (or activator concentration in the SNF) and the 
transcription rate of the repressor to determine which sets of 
parameters yielded oscillations. While the SNF, NNF and PNF 
structures have almost the same behavior with their nominal 
parameters (mean stoichiometry, amplitude and period, see 
Supplementary Figure 5A), the NNF structure oscillated over 
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the widest range of parameters and the PNF oscillated over the 
narrowest range of parameters in the simple model (Figure 5A; 
Supplementary Figure 5C). Interestingly, as the activator 
becomes more stable (i.e., the additional negative feedback 
becomes slower), the NNF structure allows sustained oscilla- 
tions over a wider range of parameters (Supplementary 
Figure 5D). Indeed, the clearance rate of the activators is 
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Figure 5 The NNF structure oscillates over the widest range of parameters. (A) The transcription rate of the repressor and the activator are changed from their 
initial value, and the range of parameters where the rhythms persist is shown. Here, dissociation constant, /C d = 10 " 5 and clearance rate of activator, 5 = 0.2. When 
/C d is varied, the NNF model still has the widest range of parameters (Supplementary Figure 5C). When 5 increases, the range of parameters that generate the sustained 
rhythms decreases (Supplementary Figure 5D). (B) Repeats the tests with the detailed mammalian model and (C) uses a Drosophila model. Details about these plots, 
as well as our methods for generating them are described in Supplementary Figure 5A-C. 



significantly slower than other circadian clock components 
(Supplementary Table 4) (Kwon et al, 2006). 

We also checked the role of the NNF structure in our detailed 
mammalian clock model. We modified the NNF structure of 
the detailed model to that of an SNF by fixing the activator 
(BMAL, CLOCK and NPAS2) concentration to the average 
value found in their WT simulations. We also constructed the 
PNF structure by converting the repressor (REV-ERBs) to an 
activator (e.g., the RORs) in the NNF structure. This did not 
significantly change the rhythms in the core feedback loop 
(Supplementary Figure 5B), matching previous studies that 
showed that the loss or change in rhythms in the activators had 
little effect on the circadian rhythms (Liu et al, 2008). It is 
tempting to conclude that the additional feedback loops 
controlling activators are not important in the circadian 
clocks. However, when we changed the transcription rate of 
the repressor {Per) and activator [Bmal, Clock and Npas2), 
the original model (with an NNF structure) had the widest 
range of parameters where oscillations occur while the PNF 
structure had the narrowest range of parameters (Figure 5B). 
Interestingly, experiments have shown that REV-ERBs play a 
more dominant role than the RORs, indicating that our 
proposed mechanism may play an important role in in vivo 



timekeeping (Liu et al, 2008). Thus, the choice of the 
additional feedback greatly affected the range of parameters 
where oscillations are seen. 

We also examined the role of the additional negative 
feedback loop in a mathematical model of the Drosophila 
circadian clock (Smolen et al, 2002). The original study that 
developed the model concluded that the NNF and SNF 
structures were equally likely to show oscillations. However, 
their study only changed transcription rates by 20%. With a 
larger perturbation of parameters, we found that the additional 
negative feedback loop significantly extends the range of 
parameters that yield oscillations (Figure 5C). 



A network design for cellular timekeeping where 
maintaining a fixed period is crucial 

The PNF structure can create a robust biological oscillator that 
has a tunable period when the additional positive feedback 
loop is fast (i.e., the activator degrades quickly) (Strieker et al, 
2008; Tsai et al, 2008; Tigges et al, 2009) (Figure 6A). 
Consistent with these findings, our simple model with the 
PNF structure has a tunable period for changes in gene 
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Figure 6 A design suitable for the cellular clocks with a fixed period. (A) A single-negative feedback loop with an additional fast positive feedback loop, with which 
activator (A) activates itself and degrades quickly. This structure has been identified in various biological oscillators like the cell cycle and pacemaker in the sino-atrial 
node (Tsai etal, 2008). (B) A single-negative feedback loop with an additional slow negative feedback loop, with which activator (A) represses itself and degrades slowly. 
Circadian clocks in mammals or Drosophila have been shown to have this structure (Supplementary Table 4) (Blau and Young, 1999; Benito etal, 2007; Liu etal, 2008). 
(C, D) The period of the NNF is nearly constant for the perturbations in transcription rates while the period of the PNF changes about two-fold. Parameters used are as in 
Supplementary Figure 5C. The period is plotted as a color where green refers to the period with the unperturbed parameters. 



expression levels (Figure 6C) . However, the simple model with 
the NNF structure has a nearly constant period in the presence 
of large changes in gene expression levels (Figure 6B and D) . 
Furthermore, this NNF structure becomes more robust as the 
additional negative feedback loop slows (i.e., the activator 
degrades more slowly) (Supplementary Figure 5D) in contrast 
to the fast positive feedback of the tunable clocks (Strieker 
et al, 2008; Tsai et al, 2008; Tigges et al, 2009). Consequently, 
our results propose two different designs for robust biological 
oscillators. The NNF structure (Figure 6B) is suitable for 
biological clocks in which the maintenance of a fixed period is 
crucial (e.g., circadian clocks). The PNF structure (Figure 6A) 
is suitable for the biological oscillators that need to tune their 
period (e.g., cell cycle or pacemaker in the sino-atrial node) 
(Tsai et al, 2008). This is also supported by mathematical 
analysis of the simple model (for more details, see 
Supplementary information) . 



Discussion 

Our work identifies several key mechanisms that allow 24-h 
rhythms in the circadian clocks of higher organisms: 
(1) Proper stoichiometric balance between the activators and 
the repressors, (2) tight binding between activators and 
repressors, (3) the NNF structure and (4) longer half-life of 
activators than repressors. These mechanisms synergistically 
generate rhythms with periods robust to gene dosages 
(Figure 6D) . The range of the stoichiometry where the rhythms 
occur widens as binding between activators and repressors 
tightens (Figure 3D). Moreover, the NNF structure regulates 
the expression of activators as well as repressors to balance 
stoichiometry (Figure 4B and C). For instance, increased 
stoichiometry (elevated repressor concentrations) strengthens 
the repression in the core negative feedback loop and reduces 
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the expression of the repressors (e.g., Pers and Crys) and Rev- 
erbs. The decreased expression of Rev-erbs weakens the 
additional negative feedback and increases the expression of 
activators (Bmall and Npas2), which lowers the stoichiometry 
(Supplementary Figure 4D). When this is done on a slower 
timescale, so that the basics of the 24-h timekeeping are 
unaffected, the robustness of the rhythms is enhanced 
(Supplementary Figure 5D). 

Relation to previous experimental data 

Many experimental observations could be interpreted as 
mechanisms by which the mammalian circadian clock 
balances stoichiometry. When the repressor (CRY) is over- 
expressed or the repressor (PER) is removed, the activator 
(BMAL1) concentration is found to increase or decrease, 
respectively (Shearman et al, 2000; Fan et al, 2007). When a 
repressor's expression is reduced, the expression of other 
repressors is increased and the expression of activators is 
decreased (Baggs et al, 2009) . Knockdown of activators yields 
opposite effects (Baggs et al, 2009). Both our detailed and 
simplified NNF models confirm these results (Figure 4B and C; 
Supplementary Figure 4D). Additionally, the rhythms of the 
mammalian circadian clock persist even after the transcription 
of all clock genes are reduced significantly (Dibner et al, 2009) . 
In agreement with these data, both the detailed and the simple 
model oscillate after significant reduction of the transcription 
rates of both activators and repressors because their stoichio- 
metry is maintained (Supplementary Figure 5E). Our study 
also suggests an underlying mechanism (ultrasensitive 
response) for a previous experimental observation showing 
that the robustness of circadian rhythms is enhanced by 
making the level of PER-CRY closer to that of CLOCK-BMAL1 
in mouse fibroblasts (Supplementary information; Supple- 
mentary Figure 3A) (Lee et al, 2011). 
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Experimental data also support the role of the slow 
additional negative feedback loop in regulating circadian 
timekeeping in higher organisms. The time course of the 
activator (BMAL1 in mammal or CLK in Drosophild) seems to 
be controlled mainly by the additional negative feedback loop 
(Rev-Erbs or Vrille) (Blau and Young, 1999; Benito et al, 2007; 
Liu et al, 2008). The elimination of additional positive 
feedback loop has little effect on circadian clocks in contrast 
to other cellular clocks based on the PNF structure (Benito 
et al, 2007; Kim and Ferrell, 2007; Liu et al, 2008; Tsai et al, 
2008). Furthermore, a key step, the clearance rate of the 
activators, which governs the timescale of the additional 
feedback loop, is significantly slower than other circadian 
clock components (Supplementary Table 4) (Kwon et al, 
2006) . Removing the slow additional negative feedback loops 
in the mammalian clock {Rev-erba ~ 1 ~ ) yields timekeeping 
where the period is not as well-maintained (Preitner et al, 
2002) . Moreover, recent studies have confirmed a pivotal role 
for the additional negative feedback loop for regulating the 
circadian rhythms via double knockout of Rev-erba and Rev- 
erbfi (Bugge et al, 2012; Cho et al, 2012). Thus, our proposed 
mechanism of robust circadian timekeeping matches known 
data on the mammalian circadian clock. Further comparison 
with known experimental data is shown in Supplementary 
Figure 7. 



Relation to previous modeling work 

Our study is the first circadian modeling study that shows the 
importance of a balanced stoichiometry in rhythm generation. 
Our results for the SNF structure match a previous model 
based on the protein sequestration (Francois and Hakim, 
2005), which focuses on other mechanisms, for example, slow 
RNA dynamics, that do not play a role in circadian clocks. We 
have identified a basic mechanism of tight binding and protein 
sequestration for generating high sensitivity, similar to what 
has been proposed in the cell cycle and synthetic studies 
(Buchler and Cross, 2009), as the key rhythm generating 
mechanism in our model. Previous circadian clock models do 
not use this mechanism, and a careful justification, based on 
experimental data from higher organisms, of the mechanisms 
for generating high sensitivity and, consequently, oscillations, 
in these models has yet to be performed (Yoo et al, 2005). In 
fact, several of these mechanisms have been called into 
question (Forger and Peskin, 2003). 

Previous models have used different mechanisms for rhythm 
generation (e.g., high-Hill coefficients) and have proposed 
different roles for the additional negative feedback loop. They 
have proposed that the additional negative feedback loop is 
capable of independent oscillations, even when the core 
negative feedback loop was removed (Leloup and Goldbeter, 
2004; Relogio et al, 2011). However, despite much experimental 
study, no oscillations have yet been found from this additional 
feedback loop in isolation (Sato et al, 2006) and the known 
phenotypes of knockout of genes in this additional feedback 
loop had not been correctly predicted (Preitner et al, 2002; 
Relogio et al, 2011). Moreover, other previous studies argued 
that the additional negative feedback loop is not important 
(Becker- Weimann et al, 2004), which does not match with 

10 Molecular Systems Biology 2012 



recent experimental data on the mammalian circadian clock 
(Bugge etal, 2012; Cho et al, 2012). We claim that the additional 
negative feedback loop is not an independent oscillator, nor 
ancillary, but acts to regulate stoichiometry. 

Interestingly, the predictions of previous modeling 
studies (Griffith, 1968; Becker- Weimann et al, 2004) match 
experimental data from the Neurospora circadian clock, in 
which a 1-1 stoichiometry is not important and the additional 
negative feedback loop seems to not play an important role 
(Baker et al, 2012). Our predictions match experimental data 
from circadian clocks in higher organisms (Supplementary 
Figures 7 and 8) . 

Proposed experiments based on model 
predictions 

Our most important prediction may be the following: when the 
stoichiometry between activators and repressors is within a 
fixed range, oscillations are sustained, and outside this range 
oscillations are damped (Figure 3). This can be tested by 
measuring the relative concentration of activators and 
repressors in many tissues and in the presence of several 
possible mutations that lead to damped or sustained rhythms. 
This has been done in WT fibroblasts and liver (Lee et al, 2001 , 
2011), but has not been done in other tissues or mutants. 
Moreover, we note that these previous experiments were done 
in population cell assays, whereas single-cell measurements 
may be needed to determine whether damped oscillations are 
the result of damped rhythms in single cell, or greater 
population desynchrony (Welsh et al, 2004; Leise et al, 2012). 

The behavior of isolated SCN neurons is similar to 
fibroblasts in that mutations of circadian genes can easily lead 
to arrhythmicity (Liu et al, 2007). We note that intercellular 
coupling in the SCN not only synchronizes SCN neurons, but 
also increases transcription of perl and per2 (Yamaguchi et al, 
2003), which may balance stoichiometry and help sustain 
rhythms when repressors are effectively removed (Tables I 
and II). Thus, we predict that increasing transcription of perl 
and/or per2 could enhance rhythmicity in isolated SCN 
neurons similar to what is seen in fibroblasts (Lee et al, 
2001). Moreover, our model predicts that cells with low 
stoichiometry (e.g., isolated SCN neurons) shows larger phase 
shifts in response to light than cells with 1-1 stoichiometry 
(e.g., SCN slices) (data not shown). It would be interesting 
future work to see whether different cell types have different 
PRCs depending on their stoichiometry. 

We also predict that tight binding between activators and 
repressors is required for rhythmicity (Figure 3D). Several 
studies have identified binding sites for PER and CRY on 
BMAL1 and CLOCK (Sato et al, 2006; Langmesser et al, 2008; 
Ye et al, 2011). Point mutations in binding sites can generate 
different binding affinities between PER-CRY and BMAL1- 
CLOCK. Comparing the experimentally measured binding 
affinities of these mutants, with the resultant rhythms, or lack 
thereof, would directly test this prediction. 

Loss of the additional negative feedback loop (e.g., in the 
Rev-erbs ~ 1 ~ , constitutive expression of Rev-erbs or constitu- 
tive expression of BMAL) is predicted to cause the intracellular 
circadian clock to oscillate over a much narrower range of 
conditions (Figure 5) . It would be interesting to test whether 
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these cells would have less temperature compensation or 
would lose rhythms more easily when other genes are knocked 
out (e.g., Cry2 ~ 1 ~ , Perl ~ 1 ~ ) . Moreover, we predict that in the 
Rev-erbs _/_ , rhythms persist in cell types with a balanced 
stoichiometry, but not in poorly balanced cells (Figure 5) . It 
would be interesting future work to investigate whether SCN 
and peripheral clocks have different phenotypes of Rev-erbs ~ f ~ 
depending on their stoichiometry. We also predict that Rev- 
erbs _/_ cells show a wider period distribution than WT 
(Figure 6) . 

Our modeling and analysis also predict that relatively stable 
activators (e.g., BMAL1 and CLOCK) in the additional negative 
feedback loop allow rhythmicity over a wide range of 
conditions (Supplementary Figure 5D). These activators can 
be destabilized with point mutations (Sahar et al, 2010) . Simply 
destabilizing the activators might lead to lower activator 
concentrations and unbalance stoichiometry, which is also 
predicted to reduce rhythmicity. However, we predict a loss of 
rhythmicity when these activators are destabilized, even when 
the overall activator concentrations are controlled for. 

Perhaps the most direct way to test our model is to build the 
clock described in our simple NNF model using the tools of 
synthetic biology. Other synthetic clocks have been built, and 
the design we propose is not more complex than what has been 
previously built (Strieker et al 2008; Tsai et al 2008; Tigges 
et al, 2009) . Validation could first be done in an analog electric 
circuit, even though this might be much less convincing. 
Building a synthetic clock would be of particular importance 
since it would be the first synthetic clock predicted to have a 
tightly regulated period. 

Future work 

Further work should explore the role of the NNF structure in the 
presence of molecular noise (Forger and Peskin, 2005; Li and Li, 
2008). Here, we studied the role of an additional negative 
feedback loop controlling the activators of the circadian clocks 
of higher organisms. Future work could consider the functions 
of the additional negative feedback loops in other organisms. In 
particular, the plant circadian clock has a different feedback 
loop structure than the mammalian or Dwsophila circadian 
clocks (Pokhilko et al 2012) . It would be interesting to see if our 
ideas carry over to other organisms and other cellular clocks. 
Furthermore, other types of feedback loops in the circadian 
clocks of higher organisms could be explored. Here, we found 
that balancing stoichiometry properly might be a universal 
principle of biological timekeeping. This finding not only is in 
agreement with experiment data from the circadian clocks in 
higher organisms, but even in agreement with the circadian 
clock in cyanobacteria as well (Rust et al 2007) . It would be 
interesting to test the role of stoichiometry in other cellular 
clocks, such as developmental clocks. 

Materials and methods 

Modifications and extensions of the detailed 
model 

The modifications and extensions of the detailed model from the 
original model (Forger and Peskin, 2003) are listed. See Supplementary 
information and Supplementary Tables 1 and 2 for details. 
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(1) Detailed modeling of additional feedback loops: The new model 
includes the secondary feedback loops, which regulate transcrip- 
tion of genes with a RORE in their promoters, including Bmals, 
Npas2 and Cryl (Preitner et al 2002; Debruyne et al 2006; Liu 
et al 2008). 

(2) Updated mechanisms of BMAL-CLOCK/NPAS2 repression: 
Matching recent findings, we updated the mechanisms by which 
the repressor (PER/CRY) inhibits the activator (BMAL-CLOCK/ 
NPAS2) (Kondratov et al, 2006; Dardente et al, 2007; Chen et al, 
2009; Yeet al, 2011). 

(3) Accounting for the heterogeneity of different genes with E-boxes: 
We introduced three different types of E-boxes for Perl /Per2/Cryl , 
Cry2 and Rev-erbs, matching experimental data (Ueda et al, 2005; 
Lee et al, 2011). 

(4) Inclusion of kinase GSK3|3: The new model includes another 
important kinase GSK3P for post-translational modification 
of the circadian clock as well as CKIe/5 (Iitaka et al, 2005; Yin 
et al, 2006). 

(5) Improved description of the effect of light on the circadian clock: 
We included a previous model of the effect of light on the circadian 
clock (Kronauer et al, 1999). 



Variables and equations of the detailed model 

The monomer proteins considered in our model are PERI, PER2, 
CRY1, CRY2, BMALs, CLOCK/NPAS2, REV-ERBs, CKI and GSK3|3. 
Although only 10 monomers are considered in the model, they can 
produce many complexes depending on the state of binding, 
phosphorylation and subcellular locations. To describe these all 
complexes, 181 variables are needed (Supplementary Tables 1 and 2): 
159 variables are for protein complexes, 12 variables are for mRNAs, 
8 variables are indicator of the promoter activity and 2 variables 
are for light effect and GSK3P activity. The reactions between 
these variables are described by ODE systems using explicit mass 
kinetics as in the original model (Forger and Peskin, 2003). More 
details of the detailed equations are provided in the Supplementary 
information. 



Parameter estimation of the detailed model 

While the original model used 36 parameters, the new model has the 
75 parameters due to the extensions and modifications of the model. 
Despite the increased number of parameters, we could get tighter 
restriction on the range of parameters with newly published data 
(listed below). Over these ranges, parameters are estimated by fitting 
to more various types of data: time courses of gene expressions and 
proteins, abundance of proteins and mutation phenotypes. 

(1) We choose 14 parameters (degradation rate of mRNAs and 
proteins) matching published experimental data. These parameter 
values were allowed to vary up to 50% from the experimentally 
determined values to account for experimental error and cellular 
heterogeneity. 

(2) PERl's phosphorylation rate is set lower than that of PER2 (Lee 
et al, 2001) . Light induced-Perl transcription is set lower than light 
induced-Per2 transcription (Challet et al, 2003). 

(3) The dissociation constant between BMALs-CLOCK and CRY is 
set greater than that between BMALs-CLOCK and PER (Chen et al, 
2009). 

(4) The ratio between cytoplasm and nucleus volume are limited to 
between 1 and 3.5 (Miller et al, 1989). 

(5) The other parameters are also restricted into a biologically 
reasonable range (see Supplementary Table S3). 

Within these restrictions, a simulated annealing method (SA, a 
global stochastic parameter searcher) (Gonzalez et al, 2007) was used 
to estimate the parameters in two steps. First, we found parameters 
that provides a good fit with mRNA and protein time profiles measured 
in mouse SCN (Reppert and Weaver, 2001; Ueda et al, 2005) and 
relative abundance of clock proteins measured in mouse liver (Lee 
et al, 2001) and fibroblast (Lee et al, 2009, 2011) (Figure 2A-C). In this 
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fitting, we used a similar cost function to that used in estimating the 
parameters of the original model (Forger and Peskin, 2003) 



10 Rj (<:..-p..\ 2 
^;=li=l U J k 

Here, j runs through 6 mRNAs and 4 proteins. Jij is the number of data 
points (12 for mRNA and 13 for protein). Sy and e t j are simulated time 
courses and experimentally measured time courses, respectively. Sy are 
normalized, in the same way as was done in the experimental data (see 
Figure 2 for details). Wy = 5 when e fJ -= 1 and Wy = l otherwise, so that 
the cost function has more weight at the peak time than other times. 
pm k and p k are maximum value of protein abundance, respectively. 
pm k and p k are normalized, so that the maximum abundance of the 
CRY1 protein is 1. 

After the first round of SA, we found several parameter sets 
qualitatively matching with experimental data on phenotypes of 
mutations of mice (WT, short, long and AR) (Table I). Then we used 
these parameter sets as initial parameter sets for another round of SA to 
get the final parameter set, which shows a quantitatively good fit with 
knockout mutation phenotype as well as time profiles (Supplementary 
Table 3). The cost function used for the second round is followed. 



£X>^-^+ £ (pm k -p k f + f£ (m Pl / mi - If + £ (ma n f 
\i=n=i lL i k V i n 

mpi and m t are simulated period and experimentally measured 
period of rhythmic phenotypes of mutations, respectively. ma n are 
simulated relative amplitude of arrhythmic phenotypes of mutation 
(e.g., Per2 _/ " or Bmall _/ ~).' 



Simulation of the detailed model 

All the simulations and parameter search were done with 150 x 8 Ghz 
CPU using MATHEMATICA 8.0 (Wolfram Research). 

Model description and analysis of the simple 
model 

The simple model is available in SBML, Mathematica, Matlab, C+ + 
and XPPAUT format from the ModelDB (Access code: 145800) (Hines 
et al, 2004; Mendes et al, 2009). The detailed model is available in 
Mathematica, Matlab and XPPAUT format from the ModelDB (Access 
code: 145801). (Schmidt and Jirstrand, 2006). See Supplementary 
information. 



Supplementary information 

Supplementary information is available at the Molecular Systems 
Biology website (www.nature.com/msb). 
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Validation of the detailed model with experimental 
data 

Time profile of mRNA and proteins 

As previously mentioned, we fit the model simulations with time 
profiles of clock mRNA and proteins in SCN to estimate the unknown 
parameters. We followed the same experimental procedures used to 
measure time profiles. The model was entrained under 12h-12h 
light/dark cycle with 500 lux light strength for 20 days. Then, the 
concentrations of mRNAs were measured during the following 48 h in 
darkness and measured time courses were compared with experi- 
mental data (Ueda et al, 2005) (Figure 2A). In the same way, the 
simulated protein time profiles are also fit with the data (Reppert and 
Weaver, 2001) (Figure 2B). 



Relative abundance of proteins 

Relative abundance among core clock proteins were compared with 
liver (Lee et al, 2001) and fibroblast data (Lee et al, 2011) because SCN 
data has not yet been reported (Figure 2C) . 



Knockout mutation phenotype 

We also tested whether our model could predict the phenotype of 
mutations of clock genes (Table I). Overall, the model simulation well 
matches with SCN or behavioral phenotype. Homozygous and 
heterozygous knockouts were simulated by reducing transcription 
rates by 100 % and 50 % , respectively. To simulate the Rev-erba ~ 1 ~ , we 
also reduced the transcription rate of the Rev-Erbs by 50%, which 
represented both Rev-erba and Rev-erbp in our model. To model the 
Bmall ~'~ , we reduced transcriptional rate of Bmals by 95%, which 
accounts for the low levels of Bmal2 when compared with Bmall (Ko 
et al, 2010). For the Clock A19/ + , the half of WT CLOCK proteins were 
mutated to be transcriptionally inactive, yet still competed with the 
remaining WT CLOCK proteins. For the CKlz tau/tau , we increase the 
CK1 phosphorylation rate for PERI and PER2 by four times. 
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